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ABSTRACT 

Motivation: Characterizing the capabilities, criticalities and response 
to perturbations of genome-scale metabolic networks is a basic 
problem with important applications. A key question concerns the 
identification of the potentially most harmful knockouts. The integra- 
tion of combinatorial methods with sampling techniques to explore the 
space of viable flux states may provide crucial insights on this issue. 
Results: We assess the replaceability of every metabolic conver- 
sion in the human red blood cell by enumerating the alternative paths 
from substrate to product, obtaining a complete map of the potential 
damage of single enzymopathies. Sampling the space of optimal flux 
states in the healthy and in the mutated cell reveals both correlations 
and complementarity between topologic and dynamical aspects. 
Availability: Related material (stoichiometric matrices, C codes, opti- 
mal flux configurations and the detailed structure of alternative paths) 
is available from http://chimera.roma1.infn.it/SYSBIO 
Contact: andrea.demartino(a)roma1 .infn.it 



1 INTRODUCTION 

Understanding metabolic activity from the underlying genotype is 
one of the most addressed problems in computational biology. Func- 
tional modularity of metabolic networks suggests that topological 
aspects may provide a key to identify a class of 'critical' structures 
or 'essential' metabolic pathways (Mahadevan and Palsson 2005; 
Samal et al 2006). However the metabolic genotype only constitu- 
tes the frame on the top of which the dynamic phenotype is built. 
Criticality of a metabolic pathway will in general depend on both the 
network reconstruction from genomic information and the 'model 
of metabolism' defined on it. In E.coli, phenotypical essentiality 
of metabolic genes has been associated with a reduced allowed 
variability of the corresponding fluxes, suggesting that dynamically 
stiff reactions may constitute an evolutionarily robust backbone of 
metabolism conserved over different species (Martelli et al. 2009). 

Here we attempt a more thorough integration of topological and 
dynamical views that will provide further insights into the efficiency, 
robustness and responsiveness of a metabolic network. We will first 
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associate the criticality of a reaction with its topological replaceabi- 
lity by enumerating the alternative paths from substrate to product 
along the network edges. Then we will validate and compare these 
results with the metabolic phenotype that results from the definition 
of a general constraint based model for metabolic flux prediction. 

We carry out our analysis on the metabolic network of 
the human red blood cell (hRBC), one of the most studied 
organisms in systems biology, from the earliest mathematical 
models of single biochemical pathways (Rapoport a/. 1976; 
Holzhiitter a/. 1985) to the currently available genome-scale 
reconstructions (Jamshidi and Palsson 2006). The reason for this 
choice lies essentially in its limited size. On one hand, it allows 
to compute reaction replaceabilities by a suitable modification of 
Johnson's exact algorithm for counting loops in a directed graph 
(Johnson 1975). On the other, it allows for the efficient applica- 
tion of various sampling methods to the space of viable flux states 
(Wiback^^a/. 2004; MarteUi a/. 2009). This is vital to address 
many important properties of erythrocytes. Indeed for some orga- 
nisms under certain conditions it is reasonable to assume that the 
metabolic activity is aimed at maximizing a subset of the metabolic 
reactions (or a function of them) associated with a certain biolo- 
gical function. In such cases the relevant flux configuration can 
be computed by standard optimization algorithms. For example, 
E.coWs metabolism has been shown to maximize biomass pro- 
duction under evolutionary pressure (Ibarra et al. 2002), but after 
a genetic knockout it responds with a minimum rearrangement of 
fluxes (Segre et al. 2002). While the production of the cof actors 
ensuring the maintenance of osmotic balance and the release of 
oxygen may be argued to be their metabolic goal, erythrocytes 
do not generically allow for such a simplification. Information- 
rich directions in flux space must be retrieved by coupling the 
underlying constraints on fluxes with other types of analyses. 
Much understanding has indeed been obtained from the uniform 
sampling of feasible states (Wiback et al. 2004; Price et al. 2004; 
Braunstein et al. 2008) and by functional studies, like the computa- 
tion of extreme pathways (Wiback and Palsson 2002), of metabolic 
regulatory structures (FricQ et al. 2003; Barretts? a/. 2006) and of 
metabolic pools (Kauffman et al. 2002). 

These aspects combined make hRBCs a key benchmark for both 
theories of metabolism and computational tools. 



© Oxford University Press 200x. 



1 



De Martina et al 




Fig. 1. Bipartite graph representation of a reaction network, with circles 
(resp. squares) denoting reactions (resp. metaboHtes). Here, reaction R uses 
metabolites a and a' as substrates to produce metabolite b. If R is knocked 
out, the conversion of a to 6 is still permitted by the alternative pathway 
a^c^d^e^b. When R is fictitiously reversed, this chain forms 
a directed loop of length 5 reactions, formed by R and by a path passing 
through £ = 4 other reactions. 

2 APPROACH 

Given a reaction network, we want to compute, for any pair of 
metabolites a and b that are respectively substrate and product 
in a reaction i (this situation will be indicated by a A 6), the 
number Af'^l^^^^t) of alternative pathways, excluding reaction z, of 
length £ allowing for the conversion a^b, see Fig.l. The idea is 
that a metabolite conversion a for which A/'^!^5(£) (or, more 
properly,^^ A/'^^j,(£)) is large will be more easily substituted in 
case of a knockout than one for which the above quantity is small. 
Finding paths connecting two points of a directed network is a 
long-studied problem in computer science. The focus is usually on 
locating the shortest paths or the fastest way to find any path. Cha- 
racterizing all the distinct paths between two vertices is however a 
less confronted issue. In our case it is crucial to avoid overcoun- 
ting, e.g. due to self-intersecting paths. Therefore we shall resort 
to an exhaustive algorithm. We will identify the substitutive paths 
using the following trick: for each pair (a, b) of metabolites such 

that a-^b, revert i fictitiously. This results in a new graph where an 

auxiliary edge 6--^a replaced the edge a-^b, see again Fig.l. Coun- 
ting the number of alternative reaction chains producing b from a 
then comes down to computing the number of directed cycles, i.e., 
non self-intersecting directed closed paths along the edges of the 

new graph, passing through the fictitious edge b --^ a. Thanks to 
the limited size of the hRBC network it is possible to solve this enu- 
meration problem exactly via Johnson's algorithm (Johnson 1975), 
briefly described in the following section. A/^!^^(£) can now be 
trivially inferred. For simplicity, £ will denote here the number of 
reactions in the alternative pathway (£ = 4 in Fig.l). 

The space of viable fluxes will be defined through a constraint- 
based approach which relies on more general assumptions than 
flux-balance analysis (FBA, (Kauffman et al. 2003)). FBA is the 
standard method to model steady-state reaction networks where 
mass balance constraints are imposed to every metabohte. For a 
reaction network with N reactions and M metabolites, let us denote 
by A and B, respectively, the M x matrices of output and input 



stoichiometric coefficients. The stoichiometric matrix is given by 
S — A — B. Letting u = {h'i)f^i denote a vector of fluxes (with 
properly chosen bounds v^^^ < < t^T^^)^ the concentrations 
c = (c")^i of metabolites vary in time according io c = Su — u, 
where u = (u"')^=i stands for the net cellular uptake of metabolite 
a {u^ > if a is a global output of metabolism, < if a is con- 
sumed by the organism, = if a is mass-balanced). Assuming a 
steady state, the concentrations are constant in time (i.e. c = 0) and 
vectors u satisfying Su = u, or 

(A-B)iy = u, (1) 

represent flux configurations ensuring that each metabolite meets its 
production or consumption constraints at fixed concentrations. As 
N is typically larger than M, the system is underdetermined and 
feasible flux states form a convex set of dimension N — rank (5') 
embedded in the AT-dimensional space of fluxes. In absence of 
a selection criterion that allows to pick one solution out of this 
set (as e.g. a maximum biomass principle), a uniform sampling 
of the solution space should be carried out. This can be achie- 
ved effectively, albeit at a considerable computational cost, by 
Monte Carlo methods (Wiback et al. 2004; Price et al. 2004) or by 
message-passing procedures (Braunstein et al. 2008). 

Here we will consider a different but related approach 
based on Von Neumann's (VN) model of reaction networks 
(Martelli et al. 2009). In the VN framework we fix the environ- 
ment through a small set of intakes on nutrients and define a 
self-consistent flux problem where the network chooses, given a 
target growth rate, how much of the nutrients to use and which 
metabolites are globally produced. Mass balance then emerges as 
a property of the solutions for some metabolites. The equations des- 
cribing the VN model have been studied by statistical mechanics 
methods in (De Martino and Marsili 2005; De Martino et al. 2007). 
For an intuitive derivation, note that the quantities Au and Bv 
represent, respectively, the total output and the total input of each 
metabolite for a given flux vector u. Then a flux vector such that 
Au > pBu, with some constant p > 0, describes a network state 
where metabolites are being produced at a rate at least equal to p, 
since for each of them the total output is at least p times the total 
input. It is simple to see that as p increases the volume of such flux 
vectors shrinks continuously (for p = every flux vector is a solu- 
tion). In particular, there exists a value of p, representing the 
maximum metabolic production rate compatible with the stoichio- 
metric constraints, above which no suitable flux vectors exist. The 
presence of conserved metabolic pools (Famili and Palsson 2005) 
implies p^ = 1 (De Martino et al. 2009), so that in metabolic 
networks optimal fluxes correspond to the solutions of 

{A - B)u > (2) 

The solutions of (2) do not coincide with those of (1) even for 
u = 0. Interestingly, a finite volume of (optimal) flux states turns 
out to satisfy the above constraints (Martelli et al. 2009). This trait 
is at odds with both the behavior of the solutions of (2) for a 
random reaction network (where a single solution survives at p^ 
(De Martino et al. 2007)) and with the optimization that is usually 
coupled to FBA (where typically a single flux state maximizes the 
objective function), and points to the robustness of metabolic phe- 
notypes. For E.coli, in particular, the solutions of (2) have been 
shown to reproduce both the large-scale organization of fluxes and 
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Fig. 2. Scheme of the hRBC metabolic network used in our analisys. 
Squares correspond to metabolites, numbers to reactions (see Table 1). 

the individual measured rates. In addition, fluxes with the smallest 
solution- to- solution fluctuations, representing the most susceptible 
parts of the network, turn out to be strongly correlated with E.coWs 
phenomenologically essential genes (Martelli et al. 2009). 

3 METHODS 

3.1 Reconstructed network 

We consider the hRBC metabolic network studied in (Price et al. 2004), a 
map of which is shown in Fig. 2; Table 1 lists reactions and the correspon- 
ding abbreviations. The network comprises three main pathways, namely 
glycolysis (reactions 1-13), the pentose phosphate (PP) pathway (14-21) 
and the adenosine metabolism, with a total of M = 39 metabolites linked by 
N = 59 reactions: 49 internal reactions (34 of which come from the splitting 
of 17 reversible processes), 3 auxiliary fluxes to maintain the osmotic equi- 
librium and the redox state of the cell (ATPase, NADHase, NADPHase) and 
7 uptake reactions to guarantee the intake of the necessary nutrients (GLU, 
ADE, ADO, INO) and of the cytosol elements (H2O, H, P^). 

3.2 Structural properties 

Structural vulnerabilities are identified by analyzing the loop structure of 
a modified metabolic reaction network, created from the original one by 
inverting -in turn- the direction of the single reaction for which we want 
to compute the replaceability, as explained in Fig.l. The fastest known exact 
algorithm (for the worst case scenario) of this cycle enumeration problem for 
a directed graph was introduced by Johnson (Johnson 1975). We shall now 
shortly describe its key ideas, referring to (Johnson 1975) for a pseudocode. 
Given a directed graph with n vertices and e edges, the algorithm is designed 
to build non-selfintersecting paths from a root vertex r to itself, loading them 
onto stacks. The main ingredients allowing for an optimal exploration of the 
graph are (a) a smart choice of the root vertex and (b) an efficient method to 
avoid duplicating cycles and repeating searches on the same portions of the 
graph. To achieve this, vertices are initially ordered in some lexicographic 
sequence, and the algorithm only selects as roots those nodes that are the 
"least" vertex (in the initial ordering) of at least one cycle. The algorithm 
described in (Tarjan 1972) garantuees to find such vertices in 0(n + e) ope- 
rations. Moreover, to avoid self-intersections, each time a node is loaded 



XT 
JN. 


Abbr. 


Chemical reaction 


1 


HK 


GLU + ATP G6P + ADP + H 


2 


PGI 


G6P ^ F6P 


3 


PFK 


F6P + ATP FDP -1- ADP + H 


4 


ALD 


T7TAT^ /"< A OT^ TAT T A T^ 

FDP ^ GA3P + DHAP 


5 


TPI 


T\TT A T* A OT* 

DHAP ^ GA3P 


6 


GAPDH 


A OTA . "VT A TA . TA 1 OTATA/~1 . "VT A TATT . TT 

GA3P+NAD+Pi ^-^►nDPG+NADH+H 


7 


PGK 


1 OTATA/^ A TATA O TA/~1 A T^TA 

13DPG + ADP ^ 3PG + ATP 


8 


DPGM 


1 OTAT~»/^ '^OTAT^/"' . TT 

13DPG 23DPG + H 


9 


DPGase 


'^OTAT^/"' TT'^/~V O TAy^ T^ 

23DPG + H20 3PG + Pj 


10 


PGM 


3PG ^ 2PG 


11 


EN 


T-v/^ TkT^T* T T/^ 

2PG ^ PEP -1- H20 


12 


PK 


T*T^T* A TAT* TT T^XT^T* A 'T^'T* 

PEP -1- ADP -1- H ^ PYR -1- ATP 


13 


LDH 


T^AT'TA TVT A TATT . TT T A A.T A TA 

PYR -1- NADH -1- H ^ LAC + NAD 


14 


G6PDH 


/""/Cn 1 XT A TAT> /CT>/~'T i XT A TAT>TT i TT 

CjoF + NADF — ^ oFCjL + NADFH + H 


15 


PGL 


/^TA/~1T TT /~\ /^TA/^/~1 TT 

6PGL + H2O ^ 6PGC + H 


16 


PDGH 


/^T^/"*/^ TVT A TAT^ TAT ^TA TVT A TAT^TT 

6PGC + NADP RL5P + NADPH + CO2 


17 


RPI 


T» T C T* T* C T* 

RL5P ^ R5P 


18 


XPI 


TAT CT^ A7"^T^ 

RL5P ^ X5P 


19 


TKI 


'\7' CT\ T^ ^T^ OTT^ /"< A O T^ 

X5P + R5P ^ S7P + GA3P 


20 


TA 


A T* d^T* T^ A T* T^/^T* 

GA3P + S7P ^ E4P + F6P 


21 


TKII 


'\.7' Cr% T"" ATX T^/"r% A OT^ 

X5P + E4P ^ F6P + GA3P 


22 


PRPPsyn 


TA ^T^ A T^T^ T^TA T^T^ A A KT\ 

R5P + ATP PRPP + AMP 


23 


PRM 


T* i T* T* ^T> 

RIP ^ R5P 


24 


HGPRT 


TATA TATA T t'S.T T T T"R /TTA ^A TA 

PRPP + HX + H2O IMP + 2 Pi 


25 


AdPRT 


TAT^ T^TA A TAT^ T T A TV TTA /A TA 

PRPP + ADE + H2O AMP + 2 Pi 


26 


PNPase 


TXT/^ T*' T T"\7' T» 1 T» 

INO -1- Pi ^ HX -1- RIP 


z / 


iiVlr ase 


TA/TP J. C\ ^ TKTO j. P j. W 

iiVlr + rl2W — ^ iiNW + r i + rl 


28 


AMPDA 


AMP + H2O IMP + NH3 


29 


AMPase 


AMP + H2O ADO + Pi + H 


30 


ADA 


ADO + H2O ^ INO + NH3 


31 


AK 


ADO + ATP ADP -h AMP 


32 


ApK 


2 ADP ^ ATP + AMP 


33 


ATPase 


ATP ADP + Pi 


34 


NADHase 


NADH NAD -h H 


35 


NADPHase 


NADPH ^ NADP + H 



Table 1. List of reactions considered in this work, including the correspon- 
ding number in the map of Fig. 2, the abbreviation and the process. The 7 
uptake fluxes, numbered 36 to 42, are as shown in Fig. 2. 



onto a stack it is also given a "blocked" status. It was proven by Johnson that 
if a vertex v stays blocked as long as every path from v to the root vertex r 
intersects the current path at a vertex other than r, the algorithm outputs all 
cycles exactly once. By sufficiently delaying the unblocking of each of these 
vertices and by keeping track of the portions of the graph that have been sear- 
ched holding the current stack, the maximum time that can elapse between 
two consecutive cycle outputs can be reduced to 0(n + e). The same holds 
for the time window before the first cycle is delivered and for the one after 
the output of the last cycle. Hence, the total time needed to list the, say, c 
cycles of the graph is 0((n + e) (c+ 1)) . In our case, each fictitious reaction 
reversal generates a new graph, so that computing the complete substitutabi- 
lity map for the hRBC requires a time of the order 0{N{n ^ e)(c + 1)). 
For practical reasons, we perform this analysis on the bipartite metabolic 
network (as in Fig.l) rather than the reduced network of Fig. 2. This implies 
that in our case n = N -\- M. 

One can in principle consider different measures of replaceability of a 

metabolic conversion b. The quantity n^^l^^^ = T.i-^a^h^^)^ which 
counts the total number of paths alternative to i from a to 6 of any length. 
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is perhaps the most obvious option, while taking into account the fact 
that longer detours are less convenient than shorter ones from an energe- 
tic viewpoint one is lead to consider e.g. exponentially- weighted functions 

like = J2i^^Pi-^Wa^b(^)- 7^-based and W-based rankings 

of metabolic conversion are rather different. The full rankings are availa- 
ble from http://chimera.romal.infn.it/SYSBIO. Here we limit ourselves to 
identifying three key reaction groups which are independent of the replacea- 
bility measure: (a) the group of reactions such that each substrate-product 
pair involved in them can be substituted (this is putatively the part of the net- 
work that is most robust to single knockouts); (b) the group of reactions that 
cannot be substituted, corresponding to the most harmful enzymopathies; (c) 
the group of reversible reactions that are only replaceable in one direction, 
corresponding to the situation in which a conversion a ^ b can only be 
substituted in one direction in case of a knockout. Note that, for topological 
reasons, intakes are not replaceable. 

3.3 Optimal flux conflgurations 

optimal flux vectors, i.e. solutions of (2), are computed by the algorithm 
introduced in (De Martino et al. 2007) based on (Krauth and Mezard 1987). 
The idea is to modify fluxes iteratively until all inqualities in (2) are satis- 
fied. Specifically, for a fixed < p < p'^ (with p"^ = 1 in our case) 
define H(p) = A — pB and let ^"(p) denote the rows of S(p), for 
a G {!,..., M}. Let also, for each iteration step t, u{t) be the flux vector 
at step t and 

m(t) = arg min ^^(p) • iy(t) (3) 

a 

At each t, the algorithm runs as foUows. If $,^'^^\p) • v{t) < 0, update 
fluxes according to 

iy,{t + 1) = max{0, u^it) + C^'\p)} (4) 

and iterate in t. Else, if ^^(*) (p) . iy(t) > Q stop, i.e. ^{t) is a solution. 

Convergence to a solution is guaranteed for all < p < p"^ 
(De Martino et al 2007), so that p'^ can be approximated with the desired 
resolution by iterating the above process for increasing values of p. To 
ensure that solutions are well defined one can either resort to setting fixed 
upper bounds on Ui 's or, as we do, impose a linear constraint of the form 
= N on the solutions (this is equivalent to singling out one flux 
as the reference unit for the other fluxes). It is convenient to initialize the 
algorithm with a random vector Different initial points generate tra- 

jectories to different solutions at p'^ and the sampling of the solution space 
thus obtained turns out to be uniform (Martelli et al. 2009). 

As a means to characterize the shape of the solution space we employ the 
average overlap between different optimal flux vectors, defined as follows. 
Let and denote two distinct solutions of (2) and, for each flux i, let 



This quantity, called the 'overlap' between a and f3, equals 1 if flux i takes 
on the same value in solutions a and [3 and decreases as the values differ 
more and more. Averaging q^fs (i) over different pairs of solutions provides 
a measure of the allowed variability of flux i (smaller variability corresponds 
to larger average overlap), complementary to the standard deviation of the 
resulting flux distribution. The complex shape of the solution space can 
then be roughly understood by distinguishing narrower directions with lar- 
ger overlap or less variable fluxes from broader ones. It is reasonable to think 
that a cell will be more sensitive to perturbations (e.g. knockouts) of fluxes 
with larger overlap. Hence analyzing the susceptibility of the solution space 
to perturbations along the directions identified by different fluxes allows to 
extract a list of the potentially more deleterious perturbations, in analogy 
with previous work on E.coli (Martelli et al. 2009). 

3.4 Response to enzymopathies 

In order to test the hRBC network against enzymopathies, we can focus on 
two types of perturbations. One can first employ a structural criterion: the 




Fig. 3. hRBC's reaction replaceability map. Black arrows denote unre- 
placeable reactions (group (b) above); thick white-filled arrows denote 
reversible reactions that can be replaced only in the direction indicated by 
the arrow (group (c)); all other reactions are fully replaceable (group (a)). 



knockout of a metabolic conversion a ^ b that is less "substituted" is more 
likely to be deleterious for the cell than the knockout of a highly replacea- 
ble conversion. The second criterion is based on fluxes: fluxes with smaller 
allowed variability (i.e. larger overlap) in the healthy cell are more likely to 
be critical links of the network than fluxes whose value can be changed over 
a larger range without losing optimality. As is to be expected, the knock- 
outs rankings produced in these ways have a large degree of similarity, as 
reactions in the group (b) discussed above coincide with the physiologically 
most critical part of hRBC's metabolism. The simplest way to simulate an 
enzymopathy on flux i is to constrain its value below a certain upper bound 
Ui. Deficiencies can be partial, i.e., of a smaller degree, the closer Ui is to 
the upper limit of the allowed range in the healthy cell, or total if 17^ = Q. 
Such constraints cause in principle a modification of the solution space along 
the direction i which in turn cascades on the entire volume, modifying the 
optimal states of the metabolic network. 

4 RESULTS 

4.1 Structural properties 

The substitutability map derived from the loop analysis is displayed 
in Fig. 3. (For the sake of simplicity we exclude the highly repla- 
ceable currency exchange fluxes from this discussion.) The most 
replaceable core of the network lies in the PP pathway (reacti- 
ons 17-21), which constitutes the main source of NADPH, a key 
metabolite in erythrocytes that limits the accumulation of peroxides 
protecting the cell from hemolysis. The high reliability coming with 
replaceability partly explains the reason why this group of reactions 
plays a central role not just as an auxiliary pathway for glycoly- 
sis, see the following analysis of fluxes. Unreplaceable reactions are 
instead lined up along glycolysis (numbers 1,6,8-13), in the bridge 
between glycolysis and the PP-pathway (14 and 16) or in auxiliary 
modules (22, 27, 29; the ADE^AMP conversion in 25 is also not 
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replaceable being directly linked to the ADE uptake). The physiolo- 
gically most deleterious knockouts (HK, PK and G6PDH) all belong 
to this group. For instance, deficiency in the level of G6PDH is the 
basis of different types of hemolytic anemias, including favism, and 
is also linked to malaria resistance (Storey 2004). Finally, there is a 
group of reversible reactions that can be replaced only in one direc- 
tion (reaction 2, 4, 7, 15, 23, 26). The last three of these could still 
be replaced in case of a knockout if a proper medium is selected. For 
instance, if reaction 15 is knocked out, it could be substituted by an 
alternative chain of reactions provided 6PGC is externally supplied. 
This is instead not possible for reaction 4 and possibly 23 (depen- 
ding on the directionality of reaction 26), as a knockout in these 
cases would necessarily result in a net production of FDP and RIP. 

4.2 Optimal flux configurations 

The flux distribution corresponding to optimal states in the healthy 
and enzyme deficient hRBC are displayed reaction by reaction in 
Fig.4, obtained by sampling 10000 solutions of (2), while a picto- 
rial representation of the optimal flux states is given in Fig.5. For the 
healthy cell (black line in Fig.4 and top left panel in Fig.5) the large 
flux backbone is formed by the second part of the glycolysis (cru- 
cial for ATP, NADH and 23DPG production) and the PP pathway 
(NADPH production). The latter gives a substantial contribution 
to the former, not just as salvage way. The adenosine metabolism 
shows instead lower flux values. In addition to GLU, which is the 
fundamental substrate for hRBCs, the INO uptake plays an import- 
ant role as an alternative way to the PP pathway. It is worth stressing 
that these solutions imply a net production of 23DPG, the crucial 
regulator for oxygen release, which is obtained without any imposed 
constraint. This picture is strongly reminiscent of the first eigenpa- 
thway obtained by extreme pathways analysis in (Price et ah 2003), 
though the thermodynamic constraints and production requirements 
used in (Price et al. 2003), including one on 23DPG, are more strict 
than the self-consistent analysis presented here. In Fig. 6 we report 
the overlap map of the hRBC. Comparing this with Fig. 3 one sees 
that the large overlap backbone coincides to a large degree with the 
structurally most vulnerable parts of the network. Note that the over- 
lap of reactions 2, 4, 15 and 26 is larger in the direction that cannot 
be replaced, further pointing to a higher susceptibility, and that cur- 
rency reactions (31-35) belong to the most constrained part of the 
network. Revealingly, however, topological and dynamical charac- 
terizations prove to be complementary in some cases. This is seen 
e.g. from reaction 3, which is flux-constrained but also highly repla- 
ceable, so that the knockout damage is limited even in presence of a 
small allowed dynamical range. (A similar picture holds for reaction 
23). To conclude, we remind that in our framework uptake fluxes are 
optimized variables not fixed by boundary conditions. In the optimal 
state five of the uptakes have a limited allowed variability, implying 
rather severe constraints on the cell's environment. 

4.3 Response to enzymopathies 

We have simulated the most studied enzymopathies by constrai- 
ning the flux of the corresponding reaction. Generically spea- 
king, the hRBC metabolism displays a large resilience against 
partial perturbations. Indeed, we have observed appreciable diffe- 
rences in relevant cellular functions compared to the non-deficient 
case only under full enzyme deficiences, as also observed in 
(Jamshidi and Pals son 2006) within a standard FBA optimization 




Fig. 6. Large overlap backbone of hRBC: black arrows mark reactions with 
an overlap larger than 0.9 in an optimal healthy cell. 



approach. Even under the most serious enzyme deficiencies the net- 
work appears to be able to maintain the production of ATP, NADH 
and NADPH almost constant, see also (Tekir et al. 2006). We focus 
here on PK and G6PDH deletions. As shown in Fig.4, the alterations 
in the flux distributions are not particularly striking and indeed we 
do not observe global flux rearrangements on the network's scale. 
The G6PDH enzymopathy appears to only cause local changes, 
confirming the structural predictions, the overlap calculations and 
also in agreement with clinical observations (Jacobasch 2000). The 
response to PK knockout is instead more marked. The synoptic ana- 
lysis of Fig.5 shows that in general the response to the perturbation 
consisted in a drop of the GLU uptake, and in a reduction of the gly- 
colitic flux, while the Rapoport-Leubering shunt (reactions 8-9) for 
the production of 23DPG remains particularly stable, as does the 
adenosine metabolism. For the glycolytic deficiences PK and HK 
we further observe an increase of the INO uptake to sustain the PP 
pathway and allow for the second part of glycolysis, and with it the 
production of ATP and NADH, to take place. Configurations corre- 
sponding to the next most severe enzymopathies (HK, EN, PGK and 
PGM) are available from http://chimera.romal.infn.it. 

5 FINAL REMARKS 

In this work we compare two robustness measures for biochemical 
networks, one based on structural properties (the reaction replacea- 
bility), the other based on dynamical capabilities (the overlaps). The 
latter depends on both the network topology and the model defined 
on it. Within VN's frame, we found that unreplaceable reactions 
largely correspond to processes with a smaller allowed flux varia- 
bility. In such directions, reaction knockouts as well as constraints 
on the fluxes are mostly harmful. Reactions with limited (but non 
zero) replaceability tend to have instead smaller overlap, so that 
while the reaction is difficult to substitute still its flux can be lar- 
gely adjusted. In an evolutionary perspective (Vitkup et al. 2006), 
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the former pathways appear as 'frozen', and perturbations at these 
nodes will require large-scale flux rearrangement, while a mutation 
affecting the latter group may be neutral and could be preserved 
across generations. Interestingly, some reactions have both a large 
overlap and a large replaceability. These, albeit structurally robust, 
are dynamically constrained and should be considered as critical 
pathways of the metabolism as well. Integrating dynamical and 
structural characterizations thus provides a rather complete pic- 
ture of the emerging network robustness. The structural analysis 
performed here was made possible by the small size of hRBC's 
metabolic network, which has served here as a model system to 
test basic concepts and algorithms. The use of the same procedure 
on a larger network, such as E.coWs, is likely to be prevented by 
CPU time growth. However, message-passing algorithms, designed 
specifically to solve combinatorial optimization or counting pro- 
blems -albeit approximately- on graphical models, may be a suitable 
replacement (Van Kerrebroeck et al. 2008). 
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Fig. 5. Maps of fluxes in the healthy hRBC (top left) and flux rerouting in three full enzymopathies: PK (top right), G6PDH (bottom left) and HK (bottom 
right). The thick arrows are dark grey if the reaction rates diminish with respect to the healthy cell, light grey if they are higher, black if they remains the same. 
The widths of the arrows are proportional to the mean value of the corresponding distributions. 
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